#!/usr/bin/env python
"""
ising_model_powerlaw

"""


from __future__ import division
import numpy as np
import networkx as nx
import numpy.random as ra
import graph_construct as gconst
import math


def find_energy(g, S, J):
    E = 0
    for i in g.nodes_iter():
        E +=  S[i] * sum_neigh(g, S, i)
    E *= -1 * J * 0.5
    return E


def sum_neigh(g, S, i):
    k = 0
    for j in g.neighbors_iter(i):
        k += S[j]
    return k

# Assumes that the first nm belong to the first mod and so on...
def calculate_local_mag(S, N, nm):
    n = int(N/nm)
    mu = 0
    for mod in range(nm):
        mu += np.abs(np.sum(S[np.arange(mod*n, (mod+1)*n)]))
    mu = mu/nm
    return mu

# Wolff Algorithm - grow a cluster
def grow_cluster(i, cluster_spin):
    S[i] *= -1
    cluster[i] = 1
    for j in g.neighbors_iter(i):
        if (cluster[j] != 1):
            try_add(j,cluster_spin)
    return 1

def try_add(site, cluster_spin):
    if (S[site] == cluster_spin):
        if (np.random.random() < add_prob):
            grow_cluster(site, cluster_spin)
    return 1


N = 512
nm = 16
avg_deg = 14
r = 0.005
p = avg_deg/(N-1)
tot_sweeps = 10000
t_trans = 1000
realisations = 1
J = 1
temp_range = np.arange(12,15.0,0.1)
#prefix = 'modular_mods_'+`nm`+'_r_'+`r`+'_avg_k_14_N_'+`N`+'_iter_'+`tot_sweeps`
prefix = 'e-r_N_512_wolff'
sweep = 0
prop_list = ""
path = '../data/'
cluster = np.zeros(N)
S = np.ones(N)
add_prob = 0

# Generate the graph
for re in range(realisations):
    #g = nx.barabasi_albert_graph(N, edges)
    g = nx.erdos_renyi_graph(N, p)
    #p_in, p_ex = gconst.find_probability(N, nm, r, p)
    #g = gconst.construct_modular_graph(N, nm, p_in, p_ex)
    print re
    print g.number_of_edges()/N
    for T in temp_range:
        accept = 0
        add_prob = 1 - math.exp( - 2 * J / T)
        mu_time_series = []
        M_time_series = []
        i = 0
        S = np.sign((ra.uniform(0.0,1.0,N) - 0.5))
        EJ = find_energy(g, S, J)
        E_avg = 0
        sweep = 0
        while( sweep < tot_sweeps ):
            for counter in g.nodes_iter():
                cluster[:] = 0
                i = np.random.randint(0,N)
                grow_cluster(i, S[i])
            if (sweep > t_trans):
                M_time_series.append(np.sum(S))
                mu = calculate_local_mag(S, N, nm)
                mu_time_series.append(mu)
            sweep += 1
        M_time_series = np.array(M_time_series)
        M = np.mean(M_time_series)
        mu = np.mean(mu_time_series)
        M2 = np.mean(M_time_series**2)
        X = (M2 - M**2)/T**2
        prop_list += "%.2f, %.3f, %.3f, %.3f, %.3f\n"%(T, M, mu, M2, X)
        print prop_list

    fname = path + 'prop_' + prefix + '_realz_'+`re`+'.data' 
    f = open(fname,"w")
    f.write(prop_list)
    f.close()





